# We assume here that the model results have been interpolated on a
# structured grid and saved in a mat file.
import scipy.io
data = scipy.io.loadmat('data.mat')
X = data['X']
Y = data['Y']
Z = data['Z']

import numpy as np
import pylab as pl
pl.figure()
# this is useful to get a tight bounding box
pl.axes([0.11,0.05,0.875,0.9])
font = {'family' : 'normal',
        'weight' : 'normal',
        'size'   : 18}
pl.rc('font', **font)

# Display the mountain
x = np.linspace(X.min(),X.max(),1000)
h_m = 2000.0
x0  = 100.0e3
a   = 10.0e3
pl.fill(x, h_m/(1.0+((x-x0)/a)**2) , 'k')

# Contour plot of the potential temperature
theta = np.reshape(Z[6,:],np.shape(X)).transpose()
print('theta min = %f'%theta.min())
print('theta max = %f'%theta.max())
levels = np.linspace(280.0,700.0,60)
c_theta = pl.contour(X,Y,theta,levels,linewidths=1.5,cmap=pl.cm.RdBu_r)

# Now the inflow wind field
step = 2
xpos = 1
u = np.reshape(Z[2,:],np.shape(X)).transpose()
w = np.reshape(Z[3,:],np.shape(X)).transpose()
q_wind = pl.quiver(X[::step,xpos],Y[::step,xpos], \
                   u[::step,xpos],w[::step,xpos], \
		   zorder=2,color='#0B3B0B',width=0.004)

# Finally, the regions where the flow is convectively unstable
Ri = np.reshape(3.0*Z[8,:]-Z[11,:]**2,np.shape(X)).transpose()
c_Ri = pl.contourf(X,Y,Ri,(-10.0,0.0),zorder=-1)

pl.ylim((-100.0,25.0e3))
pl.show()

